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The physics of the spin glass (SG) state, with magnetic moments (spins) frozen in random orien- 
tations, is one of the most intriguing problems in condensed matter physics. While most theoretical 
studies have focused on the Edwards- Anderson model of Ising spins with only discrete 'up/down' 
directions, such Ising systems are experimentally rare. LiHoa;Yi_2;F4, where the Ho"'"^ moments are 
well described by Ising spins, is an almost ideal Ising SG material. In LiHoi:Yi_a;F4, the Ho^"^ mo- 
ments interact predominantly via the inherently frustrated magnetostatic dipole-dipole interactions. 
The random frustration causing the SG behavior originates from the random substitution of dipole- 
coupled Ho^"*" by non-magnetic Y^+. In this paper, we provide compelling evidence from extensive 
computer simulations that a SG transition at nonzero temperature occurs in a realistic microscopic 
model of LiHoa;Yi_a;F4, hence resolving the long-standing, and still ongoing, controversy about the 
existence of a SG transition in disordered dipolar Ising systems. 

PACS numbers: 75.40.Cx,75.40.Mg,75.50.Lk 



The early 1970s discovery of materials failing to de- 
velop conventional long-range magnetic order down to 
zero temperature, but displaying a cusp in the magnetic 
susceptibility signaling a transition to a state of ran- 
domly frozen spins [l[, spurred thirty years of immense 
theoretical effort aimed at understanding these fascinat- 
ing spin glass (SG) systems 0, H, [l]- In that context, 
the Edwards-Anderson (EA) model of spins interacting 
via exchange interactions Jy, which can be either ferro- 
magnetic or antiferromagnetic and chosen from a frozen 
(quenched) probability distribution function, P{Jij), has 
been the subject of innumerable theoretical studies. Be- 
cause of the added simplicity of considering Ising spins 
with only two discrete 'up/down' orientations, the EA 
Ising model has attracted particular attention. However, 
because Ising magnetic materials are quite uncommon, 
most experimental studies have targeted systems where 
the moments are described instead by isotropic (Heisen- 
berg) spins that can point in any direction [3, [1, H ■ Foi' 
an Ising description to apply, the single-ion anisotropy 
energy scale must be much larger than the spin-spin in- 
teractions. This often occurs in materials where the mag- 
netic species consist of 4f rare-earth elements such as Tb, 
Ho or Dy. From that perspective, the LiHo2:Yi_a;F4 in- 
sulator has long proven to be a remarkable material to 
explore collective magnetic phenomena H, B, B [ p- fioll • 
including SG behavior within an Ising setting [111 . Il2| . 

Because of the compactness of the spin-carrying 4f or- 
bitals, magnetic exchange and superexchange between 
Ho'^^ ions is small in LiHoa;Yi_a;F4 and magnetostatic 
dipolar interactions are the predominant Ho'^^'— Ho'^^ 
couplings. Also, since the single-ion crystal field 
anisotropy is large compared to the magnetic interac- 
tions, the Ho'^^ magnetic moments can be mapped onto 
effective Ising spins that can only point parallel or an- 



tiparallel to the c-axis of the tetragonal crystalline struc- 
ture of LiHoa;Yi_a;F4 Ignoring the small nearest- 
neighbor exchange, which does not qualitatively affect 
the physics at small x, LiHo3;Yi_a;F4 can thus be de- 
scribed by a model of classical Ising spins coupled by 
long-range dipolar interactions whose Hamiltonian is: 



H = 



cr.cr., 



(1) 



Here £> > is the scale of the dipolar interactions and 
rij = \rij\, where rij = ri — rj, with Vi and rj the posi- 
tions of Ho'^+ ions i and j. zij = rij ■ z with z parallel to 
the c-axis. ei = 1 if is occupied by a magnetic Ho'^^ 
ion and = otherwise. The Ising variable ai — ±1 
for a Ho'^''' moment pointing along ±z, respectively. De- 
pending on the relative positions of two interacting mo- 
ments, the pairwise = D{rij — 3zfj)/r^j interaction 
can be either negative (ferromagnetic) or positive (anti- 
ferromagnetic). Despite the resulting geometrical frus- 
tration, pure LiHoF4 exhibits long-range dipolar ferro- 
magnetic order below a critical temperature of Tc « 1.53 
K [nl, m, [13, Ei, Ell. As Ho3+ is progressively substi- 
tuted by non-magnetic Y'^+j Tc decreases, while random 
frustration concomitantly builds up until, for Xc ~ 25%, 
dipolar Ising ferromagnetism disappears |9l. Il5|. 

It had long been thought that a dipolar Ising SG 
state exists in LiHo2;Yi_2.F4 for x = 16.5% [9j]while for 
X = 4.5%, a mysterious antiglass state occurs 0,0 1 Per- 
haps due to quantum effects It has however recently 
been suggested, on the basis of an analysis of the nonlin- 
ear magnetic susceptibility, that a SG phase might not 
actually be reahzed in LiHoa;Yi_a;F4 for x = 16.5% [l6l |. 
Even more recent work disputes this claim fvi\ , not with- 
out having generated a debate 18, l3| • To compound this 



2 



controversy, all recent numerical studies of diluted dipo- 
lar Ising models fail to find a SG transition at nonzero 
temperature 
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20|. This is in sharp contrast with the 
long-standing theoretical expectations that a transition 
should occur in this system, just as it does in the three- 
dimensional (3D) nearest-neig hbor EA model [HI, [H, [H , 
and down asymptotically to x = 0+ [24|. The field 
is thus faced with a multifaceted conundrum: is there 
a SG phase in diluted dipolar Ising materials such as 
LiHo^Yi_^F4 SHEHEi? If not, is the SG phase 
in LiHoa;Yi_a;F4 destroyed by subsidiary interactions re- 
sponsible for quantum mechanical effects that may in- 
duce an exotic (e.g. antiglass) quantum disordered state 
01? Or, is the expectation 2J, 3| that random classical 
dipolar Ising systems ought to exhibit a SG transition. 



just as it does in 3D Ising EA model [2l|, |22, l23[ , simply 



wrong? These are important questions that pertain to 
our global understanding of randomly frustrated systems 
beyond the celebrated EA model. Here, we bring new 
light on these questions by investigating model (1) via 
extensive computer simulations. 

We used Monte Carlo simulations to study Eq. [l]for a 
lattice model of LiIIOj;Yi_a;F4. We considered a tetrag- 
onal unit (C|,j(/4i/a) space group) with lattice parame- 
ters a = 5 = 5.175 A, c 10.75 A, and with four 110^+ 
ions per unit cell located at (0,0, c/2), (0, a/2, 3c/4), 
(a/2, a/2, 0) and (a/2, 0^/4). The dipolar coupling D/a^ 
was set to 0.214 K [l5[. System sizes La x La x Lc 
with L = 6, 8, 10 and an average number A'^ of spins 
A^ = AxL'^ spins were investigated via finite-size scaling 
analysis. The dipolar lattice sum in (1) was performed by 
summing an infinite array of image spins via the Ewald 
method without a demagnetization term [26j . 

Single spin-flip Monte Carlo simulations using the 
standard Metropolis algorithm is im plem ented within a 



parallel thermal tempering scheme [27|, |28[ which has 
been shown to be highly efficient in speeding up equilibra- 
tion in glassy systems. Nt replicas at different temper- 
atures were simulated in parallel with consecutive tem- 
peratures scaled by a factor a. The temperatures ex- 
plored for each replica is T^"^ = T„iina" where Tmin was 
the lowest temperature considered and n € [0, A'r — 1], 
thus the highest temperature Tmax = TminO;^^"^ and 
Q, _ "t -.1/ Tniax / Tinin ■ The acccptancc ratio for parallel 
tempering swapping is maintained above 50%. At least 
2 X 10^ Monte Carlo steps (A'mcs) per spin were per- 
formed and the last 10^ of them were used for collecting 
statistics. More than one thousand realizations of dis- 
order (A"sampio) were considered to perform the disorder 
average. Table |T] lists the parameters used in the Monte 
Carlo simulations. 

One way to monitor the freezing into a SG state is 
to calculate the overlap q{k) of two replicas with the 
same random realization of site occupancy, with q{k) = 

Ti T.i=i.2,...,N ^f^'^T^ ex.p{ik-ri) and where crf^ and cr|^'' 
are the spins of the two replicas. A standard procedure 
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TABLE I: Parameters of the Monte Carlo simulations. 

to expose a putative SG phase transition is to consider 
the dimcnsionless (scale-invariant at the critical point) 

Binder ratio [iil,!!!,!!!,^, B = i (^3 - j^^i^) , where 

(. . .) and [. . .] denote thermal average and average over 



the K 



samples 



realizations of random dilution, respectively. 
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FIG. 1; Binder ratios as a function of temperature for x = 
6.25%(left) and x = 12.5%(right), where B is dimensionless 
and T is in K units. 

Figure [T] shows B vs temperature T for a; = 6.25% 
and 12.5%. While B for different system sizes ap- 
pear to eventually merge below a certain temperature, 
no clear crossing supporting a phase transition can be 
identified. Similar results were recently obtained [l5| . 
suggesting that no finite-temperature SG transition oc- 
curs in model ([1]). That said, unambiguous B crossings 
cannot be resolved in many Monte Carlo simulations, 
even for the 3D EA Ising model where a SG transi- 
tion is believed to occur, though, in principle, crossing 
in B should be resolved when the system size is suffi- 
ciently large 21, 22, 23, [2§]. Hence, it is perhaps pre- 
mature to conclude on the basis of results such as in 
Fig. [T] that a SG transition does not occur in model 
([T]). Interestingly, Monte Carlo studies of the 3D EA 
Ising model have found that the SG correlation length, 
S,l/L, is a more suitable scale-invariant parameter to ex- 
pose a possible finite-temperature spin freezing transi- 
tion [U, [l^, If a transition occurs, ^l/L vs tem- 
perature for different L should cross at Tgg. S^l is ex- 
pected to behave asymptotically for finite L as £,l/L = 
F[(T — rsg)^"'^^'^], where F is a universal scaling function. 
The correlation length above the freezing tempera- 
ture can be approximately determined from the Fourier 
transform of the SG susceptibility, Xsg{k) = N[q'^{k)]. 
Assuming that Xsg(fe) follows an Ornstein-Zernikc form 
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above the SG transition temperature Tgg |2ll . 123 . |23|, 
Xsg(fc) (X + S,l/L can be determined via 

= (Xsg(0)/Xsg('s) ~ "^ith chosen as the 

smallest wave vector for the finite-size system, given by 
k = l-EzjicV). A suitable form for periodic boundary 
conditions is ^ = (xsg(0)/xsg(fc) - l)i/V(2 sin(|fc|/2)), 
which we use in the following calculations. 
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FIG. 2: Correlation lengths as a function of temperature for 
X = 6.25%(top) and x — 12.5%(bottom), where ^l/L is di- 
mensionless and T is in K units. The insets present the regions 
close to the crossing temperatures. 

Figuredlshows ^/L vs T for x = 6.25% and 12.5%. A 
unique and well defined crossing is observed for both con- 
centrations, providing compelling evidence that a ther- 
modynamic SG transition at Tgg > occurs in model 

Because of the small systems we need to consider be- 
cause of computational constraints, we devised an ex- 
tended scaling scheme (ESS) appropriate for the non- 
zero mean, [Jij], of the dipolar couplings Jy to analyze 
UL,T)/L, = ^^[(1 - ns/T)iTLy/-] (3^. This 

ESS is slightly different than the one used in Ref. [30| 
for the EA model with [Jij] = 0. We parameterized the 
scaling function as F{z) = Z)m=o,i,...,4 Cm(^ - zo)'^- 

After estimating Tgg from the temperature at which 
^l/L cross, the merit function. A, defined as A = 
Emc data(-^(-^)-^/?i ~ 1)^ minimized to obtain the 
coefficients c™, zq and the exponent I /v. Figure [3] shows 
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FIG. 3: Correlation lengths as a function of (1 — 
T,g/r)(rL)^''''for X = 6.25%(left) and x = 12.5%(right), 
where £,l/L is dimensionless and T is in K units. 



^l/L vs the scaling parameter z ^ {I - Tsg/T){TL)'^/'' 
where T,g = 0.047 K, 0.109 K for x = 6.25% and 12.5%, 
respectively, determined from the temperature where the 
^l/L vs r curves cross in Fig. [2] One finds the scaling ex- 
ponent l/v « 0.776 and 1/v ?» 0.782 for x = 6.25% and 
12.5%, respectively. These values are off from 1/z^ « 0.37 
for the 3D EA Ising model with [Jij] = estimated using 
an ESS with (1 - T^JT'^){TL)'^/'' as scaling parameter 
[30| . One might have expected the critical exponents of 
the dipolar model (1) to be the same as that of the 3D EA 
model, hence signaling a common universality class [2^. 
It is likely that the simulations of model ^ have not 
yet entered the asymptotic finite-size scaling regime per- 
haps. This is, in part, because of the proximity to the 
ferromagnetic phase ai x > Xc and the highly spatially 
anisotropic nature of the LiHo2;Yi_2:F4 tetragonal unit 
cell, which would both introduce corrections to scaling 
not incorporated in F{z). 
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FIG. 4: Correlation lengths estimated from Xsg(27rx/(aZ/)), 
^L,a/L, as a function of temperature for x = 6.25%(left) and 
X = 12.5%(right), where ^L,a/L is dimensionless and T is in K 
units (only data near the crossing temperatures are shown). 

We now turn to the issue of anisotropic unit cell of 
LiHo3;Yi_2;F4. The Ornstein-Zernike form for Xsg is at 
most asymptotically correct. The smallest wave vec- 
tor available in our simulation is along the c-direction 
with k = 2'Kz/{cL). However, since LiHo2;Yi_a:F4 is not 
isotropic, it is reasonable to expect that the correlation 
lengths calculated along other directions are not the same 
as that along the c-direction. Fig. 2] shows the correla- 
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further the investigation of the microscopic origin of the 
antiglass state in LiHoxYj-xFV (x = 4.5%) 0, assum- 
ing that it really exists [l9|, ISlJ ■ 
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FIG. 5: Correlation lengths and Binder ratios as a function of 
iVMcs for X = 6.25%(left) and x = 12.5%(right) with L = 10 
at the lowest temperatures. 



tion length estimated from Xsg{k = 2TTx/{aL)). We can 
clearly identify crossings at T = 0.034 K and T = 0.080 
K for X = 6.25% and x = 12.5% respectively, which 
are slightly lower than that from Xsg(fe = In z / {cL)) . 
We conjecture that, since the couplings among dipoles 
are stronger along the c-dircction then in the a-dircction 
for the LiHoi:Yi_2;F4 structure, the correlations are en- 
hanced in the former direction. This would, for small sys- 
tem sizes, move the ^l/^ crossings to a relatively higher 
temperature than for S,L,a/ L. Here too, important finite- 
size corrections are likely at play. However, without ac- 
cess to much larger system sizes and without a detailed 
analysis of the functional form of Xsg, it is impossible to 
explore this anisotropy issue further. 

The failure of some recent Monte Carlo studies [isl. [20j 
in identifying a Tgg > transition in model ^ is mainly 
because the diluted dipolar system is close to its lower 
critical dimension, as is the 3D EA model, and because of 
the sole consideration [l^ of B as an indicator of Tgg ^ 
as opposed to the more sensitive S,l/L. In addition, it is 
difficult to attain equilibrium down to the lowest temper- 
ature because of the exceedingly slow dynamics. In Fig.O 
we show the correlation lengths and Binder ratios as a 
function of A'mcs for the largest system size and the low- 
est temperature. From the figures, we believe that the 
systems are sufficiently equilibrated for extracting rea- 
sonably accurate data. 

In summary, we studied a diluted dipolar Ising model 
of LiHoa;Yi_a;F4. The spin glass (SG) correlation lengths 
show finite-size crossing as the temperature is lowered as 
well as scaling behavior, providing compelling evidence 
for a finite-temperature SG transition in model (1). It 
would be desirable to obtain data for much larger system 
sizes to improve the finite-size scaling analysis. However, 
aside from the very slow spin dynamics upon approach- 
ing Tsg, the computational effort scales as due to the 
long-range nature of the dipolar interactions, and simula- 
tions of very large system sizes will remain prohibitively 
difficult without a better algorithm. Having established 
that a SG transition occurs in the classical model ([T]), 
and for x as small as 6.25%, one may now perhaps push 
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